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ABSTRACT 

We present an analysis of the distribution of structural properties for Milky Way-mass 
halos in the Millcnnium-II Simulation (MS-II) . This simulation of structure formation 
within the standard ACDM cosmology contains thousands of Milky Way-mass halos 
and has sufficient resolution to properly resolve many subhalos per host. It thus pro- 
vides a major improvement in the statistical power available to explore the distribution 
of internal structure for halos of this mass. In addition, the MS-II contains lower res- 
olution versions of the Aquarius Project halos, allowing us to compare our results to 
simulations of six halos at a much higher resolution. We study the distributions of 
mass assembly histories, of subhalo mass functions and accretion times, and of merger 
and stripping histories for subhalos capable of impacting disks at the centers of halos. 
We show that subhalo abundances are not well-described by Poisson statistics at low 
mass, but rather are dominated by intrinsic scatter. Using the masses of subhalos at 
infall and the abundance-matching assumption, there is less than a 10% chance that 
a Milky Way halo with M V - 1T — 10 12 M Q will host two galaxies as bright as the Mag- 
ellanic Clouds. This probability rises to ~ 25% for a halo with M vir = 2.5 x 10 12 M . 
The statistics relevant for disk heating are very sensitive to the mass range that is 
considered relevant. Mergers with infall mass : redshift zero virial mass greater than 
1:30 could well impact a central galactic disk and are a near inevitability since z = 2, 
whereas only half of all halos have had a merger with infall mass : redshift zero virial 
mass greater than 1:10 over this same period. 
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1 INTRODUCTION 

Our understanding of galaxy formation is strongly influ- 
enced by our own Galaxy. Several of the apparent points 
of tension between observations and the predictions of the 
now-standard A Cold Dark Matter (ACDM) structure for- 
mation model - for example, the abundance of dark matter 
subhalos in comparison to the relative paucity of observed 
Local Group satellites and the existence of the thin Galactic 
disk in light of the ubiquity of dark matter halo mergers - 
are most clearly demonstrated from Milky Way data. While 
it is clearly an excellent testing ground for a wide array 
of phenomena related to galaxy formation, the Milky Way 
(MW) is nonetheless a single object and may not be rep- 
resentative in its properties. This is particularly true with 
respect to infrequent, stochastic phenomena such as major 
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mergers: while ACDM robustly predicts the mean number of 
mergers onto halos of a given mass, the halo-to-halo scatter 
is large. 

The situation is similar for cosmological re-simulations 
of individual MW-mass dark matter halos. In the quest to 
understand the Milky Way's formation, numericists have 
been simulating individual MW-mass dark matter halos 
at increasingly high resolution. By investing substantial 
computational resources, the most recent generation has 
reached one billion particles per halo, with particle masses 
~1000M© dDiemand et al.|2008||Springel et al.|2008||Stadel 



|et al.|[2~009[ ). Only a small number of such simulations can 
be performed and there is no guarantee that the simulated 
halos will be "typical." 

In order to explore whether various characteristics of 
the Milky Way (e.g., its satellite galaxy population, thin 
disk, or merger history) are consistent with ACDM expec- 
tations, we must understand the statistical properties pre- 
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dieted for MW-mass halos. The advent of modern galaxy 
surveys has enabled such an undertaking observationally: 
the Sloan Digital Sky Survey (SDSS; |York et al.pOOO) has 
provided a huge sample of galaxies with stellar masses sim- 
ilar to that of the MW, although only mean halo properties 
averaged over large numbers of similar galaxies can readily 
be measured dMandelbaum et al.||2006| iMore et al.||2009|. 



Predicting the formation histories and redshift zero proper- 
ties of the halos of a representative population of MW-mass 
galaxies requires numerical simulations that combine large 
volume with high mass resolution. Resolving the internal 
structure of ~ ten subhalos per host halo requires approxi- 
mately one hundred thousand particles per host, which is not 
trivial when combined with the goal of obtaining a statistical 
sample of MW-mass halos. This would, for example, require 
10 12 particles in the volume probed by the Millennium Sim- 
ulation (MS; |Springel et al.|[2005| ), a 100-fold increase over 
the particle number actually used. 

The only computationally feasible approach at present 
is to use a smaller simulation volume combined with higher 
mass resolution. In this paper, we use the Millennium-II 
Simulation (MS-II; [Boylan-Kolchin et al.|[2009[ ), which was 
designed to meet these requirements. It allows us to investi- 
gate statistical properties of MW-mass halos, including their 
assembly histories and their subhalo content. Throughout 
this work, we also compare results for our representative 
sample to the six halos from the Aquarius simulation series 
of |Springel et al.| ( |2008[ ) . These cosmological resimulations 
of Milky Way-mass dark matter halos have approximately 
one thousand times better mass resolution than the MS-II 



(see Section 2.2 1. 

The basic properties of the MS-II and of the sample 
of MW-mass dark matter halos we investigate are given in 
Section [2] together with a convergence study for subhalos. 
Section [3] discusses the assembly history of MW-mass ha- 
los in terms both of their virial mass and of their central 
gravitational potential, and relates these to their z — con- 
centrations and angular momenta. Section[4]investigates the 
subhalo population of MW-mass halos at z = 0, including 
the subhalo mass function and whether its scatter is Pois- 
sonian, the distribution of masses for the most massive sub- 
halo in each host, and typical accretion times for subhalos 
as a function of their current mass. In Section [5] we explore 
the merger histories of MW-mass halos, focusing on mergers 
that are likely to result in galaxy-galaxy mergers. Section [6] 
contains a discussion of our results and their implications, 
while Section [7] gives our conclusions. Some properties of 
satellite subhalos are quantified in the Appendix. Through- 
out this paper, all logarithms without a specified base are 
natural logarithms. 



2 NUMERICAL SIMULATIONS 

2.1 The Millennium-II Simulation 

The Millennium-II Simulation is a very large cosmological 
Af-body simulation of structure formation in the ACDM cos- 
mology. The simulation is described fully in Boylan-Kolchin 
( 2009 1. For completeness, we review some of its salient 



et al 



mass ra p = 6.885 x 10 6 h" 1 M©, from redshift z = 127 to 
z = in a periodic cube with side length L — 100 h~ 1 Mpc. 
The Plummer-equivalent force softening of the simulation is 
e = 1 h~ l kpc and was kept constant in comoving units for 
the duration of the simulation. 

The cosmological parameters used in the MS-II are 
identical to those adopted for the Millennium Simulation 
and the Aquarius simulations: 



= 1.0, fi m =0.25, fib = 0.045, «a = 0.75, 



h = 0.73, cr g = 0.9, n s = 1, 



(1) 



where h is the Hubble constant at redshift zero in units of 
lOOkms -1 Mpc -1 , as is the rms amplitude of linear mass 
fluctuations in 8 ft -1 Mpc spheres at z = 0, and n s is the 
spectral index of the primordial power spectrum. While 
the values of these parameters were originally chosen (for 
the MS) for consistency with results from the combination 
of WMAP 1-year and 2dF Galaxy R edshift Survey data 
( Spergel et al.|2003 Colless et al.|20 01), some are now only 
marginally consistent with more recent analyses. In partic- 
ular, <7g = 0.812 and n s = 0.96 are now the preferred val- 
ues based on WMAP 5-year results combined with baryon 



acoustic oscillation and Type la supernova data (Komatsu 
et al.|p009 l. For the most part, these differences are unim- 
portant for the results presented in this paper: the abun- 
dance of 10 12 h' 1 Mq halos in the WMAP5 and Millennium- 
II cosmologies differs by less than 7% at z = 0, for example. 
A discussion about possible effects of varying the assumed 
cosmology is contained in Sec. |6.3| 

In addition to the raw simulation output of six phase- 
space coordinates and one identification number per parti- 
cle, stored at 68 outputs spaced according to equation 2 of 
Boylan-Kolchin et al.| ( |2009| ), additional information about 
dark matter structures and merging histories in the MS-II 
was computed later and savedQ 

Friends-of- Friends (FOF) groups were identified on-the- 
fly during the simulation using a linking length of b = 0.2. 
All groups at each time with at least 20 particles were stored 
in halo catalogs. FOF groups were searched for bound sub- 
structure using the SUBFIND algorithm ( Springcl e t al.|200lj ). 
All substructures containing at least 20 bound particles were 
deemed to be physical subhalos and were stored in subhalo 
catalogs at each snapshot. Each subhalo in the MS-II was 
assigned a mass M su b equal to the sum of its constituent par- 
ticles; this is the mass for subhalos we use throughout this 
paper. All subhalos also have a maximum circular velocity 



V m 



V 1 

* III 



defined via 



GM sub (< r) 



(2) 



features here. 

The MS-II follows 2160 3 



10 billion particles, each of 



The main component of a FOF halo is generally a domi- 
nant subhalo^ whose particles make up most of the bound 
component of the FOF group. FOF halos can also contain 
non-dominant subhalos ("satellites"); these typically contain 



1 Dark matter halo catalogs (both FOF and SUBFIND) and merger 

trees from the MS-II are publicly available at 

http: / /www.mpa-garching.mpg.de/galform/millcnnium-II 

This is sometimes referred to as the "central" subhalo and cor- 
responds closely to the standard picture of a smooth, virializcd 
dark matter halo 
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in total only a small fraction of the mass of the dominant 
subhalo, even in the highest resolution simulations currently 
possible (|Gao et al.|2004||Springel et al.|2008||Diemand et al] 



2008| |Stadel et al.||2009| ). 

Each FOF halo containing at least one subhalo has 
an associated virial mass and virial radius, defined as the 
mass and radius of a sphere containing an average density 
A times the critical value p c (z) Throughout this paper, we 
use A = Avir, which is derived from the spherical top-hat 



collapse model (e.g., Gunn fc Gott||l972 Bryan & Norman 
19981; we denote the virial mass M v j r and the virial radius 
R v i T . Other common choices include A = 200 (correspond- 
ing to M200 and -R200) and A = 200tt m (z) [corresponding to 
M2oom and i?200m]- At z = 0, A v i r ~ 94 for the cosmology in 
equation and M200 < M v ir < M200m for any given halo. 
For dominant subhalos, M su b is comparable to, but usually 
slightly less than, M v i r . 

Merger trees were constructed at the subhalo level by 
requiring each subhalo to have at most one unique descen- 
dant. Descendants were determined by considering all sub- 
halos at snapshot Sjv+i containing at least one particle from 
a given subhalo at snapshot Sn and computing a weighted 
score based on the particles' binding energies from Sn- The 
descendant-finding algorithm also checks possible descen- 
dants at Sn+2 to account for cases where a subhalo tem- 
porarily disappears from the catalogs, usually because it 
passes very near another, more massive subhalo. Note that 
the merger trees are not constructed to explicitly follow 
mergers of FOF groups but rather to follow mergers of indi- 
vidual subhalos; the merger of two subhalos is often delayed 
relative to the merger of their original host FOF groups by 
several dynamical times. Such subhalo mergers correspond 
much more closely to the intuitive notion of disjoint objects 
merging to form a monolithic descendant than do mergers 
of FOF groups. 



2.2 The Aquarius Simulations 



The Aquarius Project ( Springel et al. 2008 \ is a suite of 



cosmological resimulations of Milky Way-mass dark matter 
halos. Six halos, denoted Aq-A through Aq-F, were selected 
at z = from a cosmological simulation known as hMS 



(Gao et al. 20081. The hMS is a lower resolution version 



of the MS-II: it follows 900 particles in the same volume, 
with the same cosmological parameters, and with the same 
amplitudes and phases 3 ! for the initial power spectrum as 
the MS-II. The halos were selected randomly from the full 
set of Milky Way-mass halos in the hMS simulation subject 
only to a weak isolation criterion: the nearest halo with mass 
greater than half that of the candidate halo was required to 
be at least 1 h~ x Mpc from the candidate halo at z = 0. 
This isolation criterion is not particularly restrictive: only 
22% of the halos from the mass-selected sample in the MS-II 
described below fail to meet this requirement. Once selected, 
the halos were re-simulated at up to five different levels of 
resolution. The highest level, which used a particle mass of 



3 The amplitudes and phases are identical for all perturbations 
with wavenumber smaller than the Nyquist frequency of the hMS; 
higher wavenumber modes in the MS-II were randomly sampled 
from the same underlying power spectrum 



led to over 1.5 billion particles within i?200m at z = 0. This 
resolution was used for only one of the six halos but all six 
were resimulated at the next resolution level, for which m p 
is at most 10 4 /i _1 Mq. These are the "level 2" simulations 
that we compare to the MS-II halos in this paper. Subhalo 
catalogs and merger trees analogous to those from the MS- 
II were also created for the Aquarius simulations. Further 
details of the Aquarius Project are presented in |SpringeI] 
et al.l d2008b. 



2.3 Milky Way-mass halos 

Recent estimates of the mass of the Galaxy's dark mat 



ter halo range from 1 - 3 x 10 12 A'/© (Wilkinson & Evans 



1999; Sakamoto eFaT][2003l fBattaglia et al. 2005| |Dehnen 



White 2008). To bracket 



|et al.| |2006; 

these estimated values of M v i r ,MW and to assess any possi- 



Xue et al.J|2008| |Li 



UP 

08f 



ble trends with halo mass, we select all halos with 10 



< 



M v ir/[ft _1 M Q ] < 10 12 ' 5 from the redshift zero output of the 
MS-II. The resulting sample of 7642 total "Milky Way-mass" 
halos forms the basis of the analysis in this paper. With the 
mass resolution of the MS-II - m p = 6.885 x 10 6 h' 1 M Q 
- this selection criterion means that the MW-mass halos 
considered here are represented by 50,000-500,000 particles. 
The typical virial radius of these halos is i? v j r w 200 kpc, 
which is 200 times larger than the Plummer-equivalent force 
softening of the MS-II. 

An alternative to selecting halos by mass is instead 
to select by Knax, which is more observationally accessible. 
Such a selection is not without complications, however, as 
the presence of a disk and bulge is likely to significantly alter 
Knax in real galaxies. In the models of Klypin et al. (20021, 



for example, the adiabatic response of the dark matter halo 
to the formation of the MW causes an increase in Knax of 
between 19% and 100%. The virial mass of the halo should 
not be modified by these effects. 



2.4 Subhalos 

Halo substructure is a direct result of the hierarchical nature 
of halo formation in the ACDM cosmology. Properties of 
subhalos such as their abundance, spatial distribution, and 
internal structure can provide an important link between 
simulations and observations, since (at least some) subhalos 
are expected to be the hosts of satellite galaxies. One should 
consider only subhalos that are reliably resolved, however, as 
subhalo structure can be modified by numerical effects at low 
resolution. The upper panel of Fig. [I] shows the relationship 
between Vm ax at z = (denoted Vmax.o) and M Bub at z = 
(M su b,o) for all satellite subhalos in our MW-mass halo 
sample that lie within distance 7? v ir from the center of their 
host halos (black lines); the three lines show the 16 th , 50 th , 
and 84 th percentiles of the distribution. Also plotted is the 
same relation for all host halos in the MS-II (red dashed 
lines). The Milky Way sample agrees very well with the full 
simulation sample. 

The upper panel of Fig. [T] also shows the best power 
law V^ax^o — Af su b,o relation fitted over the range 40 < 
Vmax.o < SOkms^ 1 (blue dashed curve) of M aub ,o = 6 x 
10 10 (Knax.o/lOOkms- 1 ) 3 ' 23 h' 1 M e . The relation follows 
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Figure 1. A comparison of various V max — M su ^ relations for 
satellite subhalos with d su b < Rvir at redshift zero. Top: the 
median and ±lcr V max ,o — M S ub relation for subhalos of Milky 
Way-mass halos (black) and for subhalos of all halos (red) from 
the MS-II. The best power law fit to the median relation is plotted 
in blue. Bottom: A direct comparison of V max ,o — M sv b g for the 
six level 2 Aquarius halos (filled black circles) and for the corre- 
sponding six halos in the MS-II (open magenta squares); the best- 
fit relation from the upper panel is also plotted as a blue dashed 
line. The two data sets agree well for Vn ax > 25kms — , while 
the MS-II points lie systematically above the Aquarius points at 
lower Vmax due to mass resolution effects. The six central halos 
are also included in the bottom panel (upper right) and are offset 
from the satellite subhalo relation because they are not affected 
by tidal stripping. 



the smallest Vmax plotted [lOkms -1 ], the MS-II data start 
to differ systematically at V max ~ 25kms~ 1 . Note that the 
host halos follow a relationship with a similar slope but with 
larger mass at fixed V ma x- This is in part because central 
subhalos are not affected by tidal stripping. In addition, the 
way SUBFIND assigns particles to bound structures, start- 
ing with all particles bound to the dominant subhalo, con- 
tributes to the offset. Given the results of Fig. [T] we adopt 
150 particles (M sub = 1.03 x 10 9 ft -1 M©), corresponding to 
Vmax.o = 28.3 km s -1 , as the minimum number of particles 
for which we consider a subhalo well enough resolved to esti- 
mate V max - This corresponds to a cut where the V max — Af su b 
relation deviates by about 10% from the power law fit in the 
upper panel of Fig. [I] 

Subhalos in high resolution simulations can lose sub- 
stantial mass as they are tidally stripped and gravitation- 



ally shocked within the potential of their hosts (e.g. 


Tormen 


et al. 


1998 


Klypin et al.||1999| Ghigna et al.||2000 


Hayashi 


et al. 


2003 


Kravtsov et al.||2004b Gao et al.||2004 


Boylan- 


Kolchin et al.|2009 1. The putative galaxies within these sub- 



halos should be much more resilient to stripping and dis- 
ruption both because they sit in the densest regions of dark 
matter halos and because the condensation of baryons to 
the centers of dark matter halos deepens the gravitational 
potential ( | White fc Rees|1978[ ). 

This differential stripping, affecting dark matter more 
strongly than stellar components, means that the relation- 
ship between the stellar content and the instantaneous dark 
matter properties of a subhalo is generally very complicated. 
It is much more likely that the mass (or V ma x) a subhalo had 
when it was still an independent halo (i.e., a dominant sub- 
halo in a FOF group) can be related to its stellar content 
in a straightforward manner. Indeed, several studies have 
shown that many features of the observed galaxy distribu- 
tion can be reproduced in A-body simulations simply by 
assuming that the stellar mass of a galaxy is a monotoni- 
cally increasing function of the maximum dark matter mass 
ever attained by the subhalo that surrounds it. This "abun- 
dance matching" assumption (Vale & Ostriker 20041 can 



reproduce observed galaxy clustering and its evolution with 



redshift (Conroy et al. 2006 Brown et al. 2008 Conroy & 



Wechsler|2009HWetzel fc White|2010| |Guo et al.|20T0 l 

Since it is useful to know when a given subhalo first be- 
came a satellite and what its mass was at that time, we trace 
each subhalo back in time and find the point at which its 
bound mass (determined by SUBFIND) reached a maximum. 
We denote this redshift by z aC c and the corresponding sub- 
halo mass and maximum circular velocity as M(z acc ) = Af acc 
and Vmax(zacc) = Vice- The subscript "acc" is used because 
subhalos usually reach their maximum mass immediately 
before being accreted into a more massive object; in this 
sense, z aC c represents an accretion redshift]^] This allows us 
to study quantities such as the distribution of z acc and the 
average mass loss of satellite subhalos. Relationships among 



the power law fit closely for V ma x > 25kms _1 , while there 
is a strong break at low V max . This is due to mass reso- 
lution, as is confirmed by the lower panel of Fig. [l] This 
panel shows the same relation but now for all six of the level 
2 Aquarius halos (black circles) and for the corresponding 
six halos in the MS-II (open magenta squares). While the 
level 2 Aquarius data follow the power-law relation down to 



4 Even though the criterion defining z acc is that the subhalo 's 
bound mass is maximized over its history, we avoid the subscript 
"max" so as not to cause any confusion with the peak in the 
circular velocity curve Vmax, which is defined at all redshifts for 
which a subhalo can be identified 
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M acc , Vkcc, and V^n ax ,o f° r subhalos are presented in the Ap- 
pendix. 



3 MAIN HALOS 

In this section, we study some properties of the main halos 
themselves. Sec. |3.1| explores the assembly histories of Milky 
Way mass halos, both in terms of M V \ T and V max . In Sec. |3.2[ 
we investigate the distribution of halo concentrations. The 
focus of Sec. |3.3| is on halo spin parameters. Throughout this 
section, we compare properties of the six Aquarius halos to 
properties of the full MS-II sample of "Milky Way" halos. 



3.1 Halo Assembly 

The mass assembly histories (MAHs) of dark matter halos 
have been studied extensively via iV-body simulations and 
extended Press-Schechter ( 1974 \ theory (ILacey & Cole 1993 



1994 



Wechsler et al.|2002 



van den Bosch 2002; Zhao et al. 



2003b a; Tasit siomi et al.|2004||Cohn fc W hite 200 5||Neistein 
""Li et al.||2008||Zhao et al.||2009J [McBride et al. 



20091. As a direct result of this work, many basic proper- 



ties of MAHs have been established. For example, |Wechsler] 
et al.| ( |2002[ ) found that MAHs of individual halos follow a 
one-parameter family M(z) = Moe~ KZ . The parameter k 
is directly related to the typical formation time for halos of 
mass Mo and is an increasing function of mass. Tasitsiomi 
et al.l J20041 showed that this exponential form was unable 



to fit individual cluster-mass MAHs in some cases and that 
a more general function, M(z) — Mo (l + z) v e~ KZ , provided 



a better match on a case-by-case basis. McBride et al. ( 2009 1 



confirmed and extended this result over a wide range in halo 
masses. 

The median MAH from the full MW sample is plotted 
as the dashed line in Fig. [2] the gray shaded region marks 
the ±1 a region for the MAHs. A purely exponential form is 
a poor fit to the median MAH, while the modified form of 
Tasitsiom i et al.] with n — 0.795 and r\ = 0.435 provides a 
fit to within 10% for 2 < 5. The median MAH can be fitted 
extremely well by 



M{z) = Mo (1 + z) v exp[-«' (y/(l + z) - 1)] 



(3) 



with k' = 4.90 and r\ — 2.23: this fit matches the measured 
median MAH to within 1.6% for z < 10. The solid lines in 
Fig.[2]show the MAHs for the Aquarius halos measured from 
the MS-II^The Aquarius halos sample the full distribution 
of MAHs: halos A and C form early, halo F forms late, and 
halos D and E track the median MAH. Halo B trails the 
mean MAH at early times but catches up at z ~ 2. 

While the MAHs are normalized to unity at z — 0, the 
la scatter in M(zq)/M(z = 0) grows fairly quickly with 
redshift, reaching a factor of 2.5 at z = 2 and 4 at z = 4. 
The scatter in log(l + z) at fixed M/M(z = 0) is nearly 
constant for M/M(z = 0) < 0.6. As a result, formation 
times zf defined with respect to fixed fractions of the final 



L ' b igure 13 of |Boylan-Kolchin et al] |2009[ | shows that there is 
excellent agreement between the MAHs of the level 2 Aquarius 
simulations and the corresponding halos in the MS-II. 




Figure 2. Top: Mass assembly histories of Milky Way-mass ha- 
los. The median MAH is shown as a dashed curve, while the gray 
shaded region contains 68% of the distribution. Also plotted are 
the MAHs for the MS-II versions of the six Aquarius halos. The 
Aquarius halos provide a fairly representative sample of MAHs; 
halos A and C have quiescent merger histories, while halo F has a 
recent major merger. Bottom: V max accretion histories from the 
MS-II. As in the upper panel, the median relation is plotted as a 
dashed curve while the ± 1 a region is shown as the gray shaded 
region. The earliest-forming Aquarius halos in terms of MAHs - 
halos A and C — also assemble earliest as judged by V^ ax , with 
V^ ax changing by less than 25% from z = 5 to z = 0. Halo F, 
on the other hand, undergoes a major merger at z ss 0.4 that 



mass - i.e., M(z/) = f M(z — 0) - will all have similar 
scatter in log(l + zf), approximately 0.25, for / < 0.6. 

While characterizing the growth of Af v i r for a halo that 
reaches MW-like masses at z = is useful, it is far from a 
complete description of that halo's assembly history. Many 
studies have shown that halo growth occurs inside-out, with 
the central regions built up first and the outskirts assem- 
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bled later (e.g., 


Fukushige & Makino| 2001 


Loeb & Peebles 


2003 


Zhao et a 


. 2003b| |Diemand et al.||2007| |Cuesta et al. 


2008 


1. Furthermore, defining halo boundaries with respect 



to Ap c (z) means that a halo's virial mass will increase with 
time, even in the absence of any physical accretion, simply 
due to the re-definition of the boundary. A halo can therefore 
grow substantially in M v i r without any appreciable change 
in the mass at small radii, where the main baryonic com- 
ponent should lie. Accordingly, it is useful to find a way to 
characterize the build-up of the central regions of halos in 
addition to studying the growth of M v i r . 

To this end, we consider the growth of V^ax in the bot- 
tom panel of Fig. [2] For the sample of Milky Way- mass halos 
considered here, i? max is on average a factor of 6.1 smaller 
than i?vir, so Vm ax probes the mass distribution on a scales 
~ 6 times smaller than R v i T . Additionally, the gravitational 
potential energy per unit mass of a halo is proportional to 
V^ ax , so studying the evolution in V^ ax probes the growth of 
the dark matter halo's central potential. The dashed curve 
in the bottom panel of Fig. [2] show the median relation 
while the shaded region shows the ±1 a scatter. For redshifts 
z < 10, the median value of V^ ax (z) can be approximated 
to within 1.3% as 




V msx (z) = V maX! o (1 + 2) 



(4) 



Figure 3. The probability distribution of halo concentrations c 
(solid black line), along with the best log- normal fit (dashed line) 
with mean (logc) = 2.56 and standard deviation <ri ogc = 0.272. 
Values for the individual Aquarius halos in the MS-II are shown 
as colored vertical lines. See the text for a discussion of how the 
concentrations are computed. 



Qualitatively, the trend in V^ ax (z) is the same as in 
Mvir (z): both grow rapidly at early times and slowly at late 
times. At a more detailed level, however, the two show im- 
portant differences in growth. V^ ax - and consequently, the 
central potential - approaches its redshift zero value much 
earlier than M v i r : if we define a formation redshift Z f v ELS the 
time when V^ ax reaches half of its present value, then the 
median Zf ]V is approximately 4, in comparison with Zf « 1.2. 
This is a reminder that the central potential of a halo is set 
much earlier than its virial mass. The scatter in log(l + z) 
at fixed Vm ax is noticeably larger than at fixed M v i r , show- 
ing that the spread in formation times is larger if these are 
defined using V^ ax rather than M V1T . 

As in the upper panel of Fig. [2] the build-up of V^ ax 
for the Aquarius halos (the colored lines) spans the range 
given by the full sample of MW-mass halos. Furthermore, 
the behavior of mass growth for individual halos is matched 
by that of Vm ax - Halos A and C also assemble their central 
potentials at an early time, while F forms much later. 

A stable central potential is likely to be conducive to 
disk galaxy formation and evolution, as potential fluctu- 
ations may drive "secular" evolution and help transform 
disks into bulges. Halos A-E have nearly constant central 
potentials from z = 1 to z — (in the case of C, from 
z = 2), indicating they are possible disk galaxy hosts. Halo 
F, which experiences a late major merger, would likely host 
a spheroid-dominated galaxy. These expectations are gener- 
ally borne out by SPH simulations of the Aquarius halos by 
Scannapieco et al. (20091. Using a dark matter mass reso- 



lution 2 to 3 times better than that of the MS-II, they find 
that halos C, D, and E all have well-defined disks and that 
halos A and B also contain (smaller) disks, while halo F has 
a spheroid with no disk at all. 



3.2 Halo Concentrations 

While numerical simulations have shown that dark matter 



halos have a nearly universal density profile (Navarro et al 



1996 



1997 hereafter NFW), the radial scale of the profile 

1/3 

does not follow the virial scaling of R oc M vir . An addi- 
tional parameter is therefore required to specify the profile 
of the typical halo of a given M vlr \ a common choice is halo 
concentration c, defined as the ratio of the virial radius to 
the radius r_2 at which the logarithmic slope of the den- 
sity profile reaches —2. The scatter in c at fixed halo mass 
is related to the diversity of halo formation histories. Pre- 
vious studies have shown that average halo concentrations 
decrease weakly with halo mass - c oc M -01 - with a fairly 
large scatter at fixed mass, <r(log 10 M) « 0.1 — 0.14 (Navarro 
et al.|1997j|Jing|2000[|Bullock et al |2001b||Neto et al.|2007[ 



Gao et al.|2008| |Macci6 et al.|2008| |Zhao et al.|2009| >. 

We can estimate a concentration parameter for each 
halo in our MW-mass sample by assuming that each halo's 
density structure can be fitted with an NFW profile; we then 
have i? max = 2.16r_2 and therefore, c = 2.16 7i v i r /J? max . 
The measured i? max values for the Aquarius halos in the MS- 
II differ from those in the level 2 Aquarius simulations by less 
than 10%, showing this radius to be accurately determined 
at MS-II resolution. The probability distribution function 
for log c from our MW-mass halo sample is plotted as a solid 
black line in Fig. [3] The dashed line shows the best-fitting 
normal distribution with (logc) = 2.56 and <7i g C = 0.272. 
The dark (light) shaded region marks the ±1 a (±2 <r) range 
of this distribution. There is a slight but noticeable excess 
of halos at low concentrations compared to the best-fitting 
log-normal distribution; |Neto et al.| ( |2007[ ) have shown that 
this is due to a population of unrelaxed halos. 

Vertical colored lines mark the concentrations of the 
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Aquarius halos in the MS-II, computed as described above. 
Halos B, D, E, and F lie near the peak of the distribu- 
tion while halos A and C lie in the high concentration tail 
(though see the caveat about halo A below). This behavior 
is expected given the results of Fig. [2] the concentration of a 
halo is (loosely) related to the density of the universe when 
the halo formed. Early-forming halos such as A and C are 
more concentrated because the universe was denser when 
they formed. 

It is important to note that we did not fit the density 
profiles of these halos to a specific parametric form such as 
NFW when estimating the concentrations presented here. 
Instead, i? ma x and J? v i r were measured directly from the 



simulation data using SUBFIND. Navarro et al. (2010) have 



made a thorough study of the structure of the Aquarius ha- 
los at z = 0, concluding that their density profiles are not 
perfectly universal and therefore cannot be fully described 
by a model with two scale parameters but no shape parame- 



ter (such as NFW). In agreement with Neto et al. (20071 and 
Gao et al. (20081, they found significantly better fits with 



the three parameter model of |Einasto| ( |1965p , although the 
mass range spanned by the Aquarius halos is too small to see 
the systematic redshift dependence of the shape parameter 
that was measured in the earlier work. We emphasize that 
profile fitting typically attempts to represent the widest pos- 
sible range in radius but tends to focus on the inner regions, 
usually starting at the minimum resolvable radius. This is in 
part because inner regions of halos are more relaxed and also 
more regular than the outskirts. For example, Navarro et al. 
show that the Aquarius A-l halo is fitted extremely well by 
an Einasto profile over the range 0.1 h,- 1 <r< 25/i _1 kpc 
but that at larger radii, its density profile rises above this 
model due to substructure and recently accreted material. 
The latter is a substantial contribution to the mass and leads 
to an increase in the virial radius of the halo relative to an 
extrapolation of the profile fitted to the inner region. This is 
not the case for any of the other halos (see fig. 3 of Navarro 
et al.) In this sense, the measured concentration of the A 
halo is larger than it "should be." 



3.3 Halo Spin 

The angular momentum content of a galaxy is often assumed 
to be connected to the spin of its dark matter halo. In the 
commonly-adopted model of Mo, Mao, & White ( 1998 1, the 



specific angular momentum of a galactic disk is a fixed frac- 
tion of that of its host halo. Many more recent versions of 
the Mo, Mao, and White model (e.g., |Dutton et al.|p007l 



modify this assumption only slightly. Quantifying the angu- 
lar momentum content of dark matter halos may therefore 
help to connect dissipationless A-body simulations to the 
properties of disk galaxies. 

A convenient parametrization of a dark matter halo's 



angular momentum is the modified spin parameter A' ( Bui 
lock et al.|2001a| ), defined as 

1 7 



A' 



y/2 Vvir Rvir 



(5) 



where j is the specific angular momentum of the halo. The 
distribution of spin parameters for dark matter halos as 
a function of their mass has been studied by several au- 
thors ( |Barnes fc Efstathiou||1987| |Warren et al.||1992| |Cole 




-6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 
log A' 

Figure 4. The probability distribution for the spin parameter A' 
for Milky Way-mass halos in the MS-II (solid black curve), along 
with ±1 and 2 a regions (shaded dark and light gray, respectively). 
Dashed black line: best-fitting log-normal distribution, with mean 
(log A'} = —3.41 and standard deviation <J\ og y = 0.580. Dashed 
cyan line: fit to the distribution given in equation with Ao = 
0.0358 and x = 2.50. 



& Lacey 1996 Lemson &i Kauffmannj 1999| Bullock et al 



|2001a| |Bailin fc Steinmetz||2005| |Macci6 et al.||2007| |Bett 
|et al.|2007[|Macci6 et al.|2008|). One of t he most recent, and 
extensive, studies is that of Bett et al. (20071, who showed 



that the distribution of halo spins in the MS is well-fitted 
by the function 



PQog\')=A 



exp 



3/x' 



(6) 



In Figure [4j we show the probability distribution function 
of A' for the MS-II MW-mass halos (solid curve), with ±1 
and 2 a regions given by the dark and light shaded re- 
gions, respectively. The best-fitting log-normal distribution 
has (log A'} = —3.41 and cr log v =0.580 and is shown as a 
dashed black line. The dashed cyan line shows a fit to equa- 
tion |6| with Ao = 0.0358 and \ ~ 2-5; this provides a sig- 
nificantly better description of the data than the log-normal 
fit. 

The A' values for all of the Aquarius halos in the MS-II 
are shown as vertical tick marks in Fig. [4] While all of the 
halos lie within the 2 a region, five of the six lie below the 
median A'. The lone exception, halo F, had a recent major 
merger (see Fig. Recent major merger remnants tend to 
have spin parameters that are higher than typical (Vitvit- 
ska et al.|[2002| |Maller et aL"1|2002| >. The spin parameter of 



these halos usually decreases after the major merger due to 
further accretion, which brings in additional non-aligned an- 
gular momentum. It therefore seems that at least five, and 
possibly all six, of the Aquarius halos lie in the lower half of 
the A' distribution. 
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Figure 5. The differential subhalo mass function. Data from the 
2039 MW-mass halos in the MS-II with 10 12 < M^/h^ 1 Mq < 
10 12 ' 5 are shown as the black data points, while the best-fitting 
relation using equation (J8j is shown as a dotted black line. Results 
are plotted both for rcdshift zero subhalo masses (M au ^, lower 
points) and for masses at accretion (M acc , upper points; offset 0.5 
dex vertically for clarity). Extrapolation of the MS-II results to 
low fj. agrees well with results from the individual level 2 Aquarius 
simulations (colored lines). 



4 SUBHALOS 

The abundance of subhalos in MW-mass halos has been the 
topic of fierce debate since late in the 1990's, when cosmo- 
logical simulations were finally able to reliably resolve dark 
matter substructure within galaxy-scale halos. These stud- 
ies predicted that MW-mass halos should have hundreds of 
dark matter subhalos with V^ax > lOkms" 1 , while only 
of order ten such satellite galaxies were seen in the Local 
Group i |Klypin et aT7||1999| poore et al.||1999[ ). Since then, 
the discovery of several faint Local Group satellites in SDSS 
data (e.g.. 



Belokurov et al. 20071 has extended the range 



of observed satellite luminosities by almost four orders of 
magnitude, to ~ 10 3 Lq , roughly doubling the observational 
sample. On the other hand, the dynamic range of dark mat- 
ter simulations has increased by an additional three decades 
in mass and the number of resolved dark matter subhalos 
has increased by a factor of approximately one thousand: 
Springelet al] i ]2008| find 10 5 subhalos with V max > 2kms _1 
inside i?200m of the Aq-A-1 simulation at z = (see also 
Diemand et al.|2008 and Stadel et al.|2009 I . It is unambigu- 



ous that the observed number of luminous satellites with 



V m 



> 



10 falls far short of the expected number of dark 



matter subhalos above the same V ma x threshold for standard 
ACDM models. The resolution of this discrepancy likely lies 
in the strong dependence of galaxy formation physics on 
Af acc (see Kravtsov|2010 for a recent review) combined with 
the limited sky coverage and the surface brightness detec- 
tion threshold of the SDSS ( |Tollerud et al.|2008| [Walsh et al. 
|2009[ ). Substructure in galaxy-mass halos is therefore still 
of great interest for clarifying the properties both of dark 



Figure 6. The cumulative subhalo mass function for MW-mass 
halos having 10 12 < M^/h' 1 Mq < 10 12 5 . Results are shown 
both in terms of M su t, q (lower solid line) and M aC c (upper solid 
line). The dotted curves show equation |7| using the parameters 
determined by fitting the differential mass functions (Fig. [5j . 



matter and of galaxy formation at low mass (and low star 
formation efficiency). In this section, we explore the sub- 
halo content of the MW-mass halo sample from the MS-II 
at z = 0. 

In the analysis that follows, the dominant subhalo is ex- 
cluded in all cases. Furthermore, we consider only subhalos 
lying within i? v i r ; from their host; this provides consistency 
when comparing halos of different host mass. It is very im- 
portant to note that using a different limiting radius such as 
R200 or i?200m will strongly affect subhalo abundance statis- 
tics because the radial distribution of subhalos is heavily 
weighted toward large radii ( De Lucia et al.|2004 Gao et al. 
|2004| [Springel et al.|2008| . This choice also affects statistics 
measuring the dynamical evolution of subhalos, as subha- 
los falling into their host for the first time experience tidal 
stripping that depends on the local density and therefore, 
on the pericenters of their orbits about the halo. 



4.1 Subhalo mass function 



As shown by several authors (e.g., Gao et al.|2004 Giocoli 



et al. 2008 Angulo et al. 20091, the abundance of subha 



los with a given fraction of their host's mass varies sys- 
tematically, with more massive hosts having more subha- 
los. This reflects the later formation times of more mas- 
sive halos, which leave their subhalos less time to be tidally 
stripped. In the MS-II, the normalization of the differen- 
tial mass function increases by approximately 15% for each 
decade in host mass. In order to maximize our resolution, 
in this section we use only the 2039 MS-II host halos with 
10 12 < Af vir //i -1 Mq < 10 12 ' 5 when computing the mass 
function. 

The differential mass function of (non-dominant) sub- 
halos is shown for this sample in Fig. [5] Results are plot- 
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Figure 7. The cumulative subhalo velocity function for MW- 
mass halos having 10 12 < M v i r //i — 1 A/q < 10 125 (solid line), 
plotted in terms of v = Vmax/V^i,. host - The dotted line shows the 
analog of equation determined by fitting the differential ve- 
locity function; the best-fitting parameters are given in the figure. 



ted as black points both for the actual redshift zero masses 
(Af S ub; lower data points) and also for the accretion masses 
(Ma CC ; upper points, offset vertically by 0.5 dex for clarity). 
For comparison, subhalo mass functions obtained from the 
individual level 2 Aquarius simulations (thick colored lines) 
are also plotted in Fig. [5] 

We assume that the cumulative mass function can be 
represented by a power law with an exponential cut-off at 
high masses: 



{N(> /i)) 



exp 



Mcut 



(-) 



The differential mass function then also has a simple func- 
tional form: 



d{N) 
dlog^i 



Mcut 



{N(> „)) 



(8) 



Equation Q provides an excellent fit to the MS-II data for 
Jh = 0.01, ^ cut = 0.096, a = -0.935, and b = 1.29. This fit 
is plotted as a dotted line in Fig. [5] Although the range of 
overlap between the MS-II and Aquarius data is relatively 
small, extrapolating the fit from the MS-II data to low n 
agrees extremely well with the Aquarius data all the way to 
the Aquarius level 2 resolution limit, fi m 10~ 6 . 

We also use equation Q to fit the subhalo mass func- 
tion computed in terms of M acc . The comparison of MS and 
MS-II results in |Guo et aT| ( |2010[ > suggests that N p (z acc ) > 
1500 per subhalo is required to obtain converged results for 
M acc . We are therefore able to probe (iV(> ^t a0 c)) only for 

> 10~ 2 . This limited range makes the determination of 
the slope a using the MS-II data alone nearly impossible. 



Giocoli et al. (2008) have suggested that this slope is close 
to that of the redshift zero cumulative mass function, and 



we also fit ^i acc with equation Q holding a fixed to -0.935, 
the value obtained from fitting the differential mass func- 
tion for M su b. A good fit can be obtained with Jii — 0.038, 
Mcut = 0.225, and b — 0.75 and is plotted as the upper dotted 
line in Fig. [5] Again, the extrapolation of the MS-II fit to low 
masses agrees very well with subhalo mass functions com- 
puted directly from the level 2 Aquarius simulations (upper 
set of solid curves) . This shows that our adopted value of a 
is indeed appropriate. 

The corresponding cumulative subhalo mass functions 
are plotted in Fig. [6] for redshift zero subhalo masses (M su b, 
lower solid curve) and for subhalo masses at accretion (M acc , 
upper solid curve). The dotted lines show equation |7| with 
parameters taken directly from fits to the differential mass 
functions (i.e., we do not fit the cumulative mass functions 
independently). These curves are excellent representations 
of the data for occupation numbers {N(> /i)} > 0.03. At 
lower N, the presence of a very small number of ongoing 
major mergers results in an excess cumulative abundance 
compared to the fit for M 8u b. 

The abundance of subhalos as a function of Vni a x is a 
related quantity, and one that is often used when comparing 
simulation data to observations because V^n ax is less affected 
than M su b by the dynamical evolution of a subhalo within 
its host. Fig.[7]shows the cumulative abundance of subhalos, 
(N(> u)), in terms oiv = V ma x/Kd r ,host (solid curve). It ex- 
hibits the same behavior as {N(> /u)): a power law at low v 
with an exponential cut-off at high v. By fitting the differ- 
ential velocity function (to ensure that errors in each bin are 
uncorrelated) to the same functional form as in equation |8| 
and converting the fit to the cumulative velocity function, 
we find that the slope of {N(> u)) is -2.98 at low v and that 
{N(> v)) = 1 at v m 0.4; this fit is shown as the dotted 
line in Fig. [7] This \ow-v slope is in good agreement with re- 
sults from the individual, high resolution Aquarius and Via 
Lactea simulations. 



4.2 Scatter in the subhalo mass function 

The halo-to-halo scatter about the mean relation in equa- 
tion Q is a matter of considerable interest for several rea- 
sons. A number of groups have invested considerable com- 
putational resources into simulating individual MW-mass 
halos at extremely high resolution. State-of-the-art simula- 
tions currently use approximately one billion particles to 
resolve the mass distribution within M v i r (Diemand et al. 
[2008| |Springel et al.||2008"l |Stadel et aL||2009[ ). Such calcu- 



lations are extremely expensive, however, so only a small 
number can be performed. In order to interpret their re- 
sults, we must know how much scatter is expected among 
halos of the same mass. This scatter is also a fundamen- 
tal input parameter for halo occupation distribution (HOD) 
models (e.g.,|Ma fc Fry|2000| [Peacock fc Smith|2000||Seljak 
2000] |Scoccimarro et al |2001| |Berlind fc Weinberg||2002 \ 



We find that the scatter in the cumulative mass func- 
tion at fixed fj, can be well-modeled by postulating that the 
variance in A^(> /i) is the sum of two contributions, one 
due to Poisson fluctuations (op = (Af)) and one intrinsic 
(a? oc <iV) 2 ): 



2 2.2 

a = CTp + o"i . 



(9) 

This is demonstrated in Fig. [8] The solid curves show a /up 
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Figure 8. The ratio of the measured dispersion a in N(> fi) 
relative to that expected from a Poisson distribution with mean 
value equal to (N(> fi)). The dispersion is larger than Poissonian 
and the ratio of the two increases systematically to lower fi. The 
dashed gray line shows a model in which the variance is a linear 
sum of Poisson fluctuations and an 18% fractional intrinsic scatter 
[equation |9j], while the shaded gray region encompasses the same 
model with intrinsic scatter between 14 and 22%. 



as a function of ii for host halos of different masses, from 
Milky Way-mass hosts to halos corresponding to rich galaxy 
clusters. We have used the MS for host halos with M v i r > 
10 14 /i _1 Mq, as the volume of the MS-II does not provide 
sufficient statistics at this mass scale; for all other masses 
we use the MS-II. The dashed gray line shows the model of 
equation Q with a fractional intrinsic scatter si = a\/(N) 
of 18%, while the shaded gray region corresponds to varying 
si from 14% to 22%. 

This analysis shows that the scatter in the subhalo mass 
function is nearly Poissonian for massive subhalos (ii > 
5 x 10 -3 , corresponding to occupation numbers of (TV) < 2) 
but is broader than Poissonian at low masses (/i < 10 -3 or 
(TV) > 8). This behavior agrees with the mass functions com- 
puted from the level 2 Aquarius simulations, which probe an 
additional three decades in mass to pi ~ 10~ 6 ' 5 . The scat- 
ter in the subhalo abundance among these simulations is 
fairly small: the fractional scatter, defined as the standard 
deviation in TV(> pi) normalized by the mean of TV(> pt), 
is constant over the range 10 ~ 6 < pi < 10 -4 ' 5 at approx- 
imately 12%. This is somewhat smaller than our estimate 
from MS-II data but is consistent given the small number of 
Aquarius halos. Poisson scatter would give 4.5% at pt = 10~ 5 
and 1.5% at pi = 10 -6 . This is clearly inconsistent with the 
Aquarius results. 

Many previous studies have analyzed subhalo occupa- 
tion statistics using reduced moments {cij} of the subhalo 
mass function, defined via 

(TV(TV-1)...(TV- j + l)) 1/j 



otj 0) 



(TV) 



(j = 2,3,4,...) .(10) 



Figure 9. The reduced second moment parameter »2 versus 
fi for several ranges of M v i r . The behavior of ctiifJ.) is nearly 
independent of host mass: 02 = 1 + e, with e £3 0.02, for /1 < 
4 X 10 -3 . Despite this, a comparison with Fig. |8]shows that the 
satellite HOD is not Poisson at these mass fractions. 



of HOD modeling, have used the ctj values as measures of 
how close the scatter in subhalo count is to Poissonian: for a 
Poisson distribution, Qj = 1 for all j. All studies have found 
that «2 ~ 1 over a wide range in mass and have concluded 
that the scatter is indeed Poissonian (e.g., |Kravtsov et al 



2004a| |Zheng et al?]|2005| ). While the finding that a 2 « 1 
does indeed justify the replacement of (TV (N—l)} with (TV) 2 
when computing correlation functions in HOD modeling, it 
does not show that the variation about (TV) at fixed pi can be 
well-described by Poisson statistics: distributions that differ 
substantially from Poisson can result in mere percent-level 
variations from ol% = 1. 

Figure[9]shows 02 as a function of pi for the same ranges 
of host halo mass as Fig. [8] The behavior is similar at all 
masses: ct2 is very close to, but slightly greater than, unity 
for pt < 4 x 10~ 3 , corresponding to occupation numbers 
(N) > 2. The deviations from CV2 = 1 are quite small - 
2 to 3% - over this range. This is precisely the range where 
the scatter in 7V(> y,) deviates systematically and substan- 
tially from the Poisson expectation: Fig. [5] shows that the 
scatter is 30% broader than Poissonian at n = 4 X 10~ 4 and 
nearly twice as broad at fi — 10~ 4 , whereas = 1.02—1.025 
over this same mass range. The satellite HOD is certainly 
not Poisson for fi < 10 -3 even though 012 is within 2% of 
unity. While Poisson scatter is an excellent approximation 
for (N) < 3, it is increasingly inaccurate at higher occupa- 
tion number. 

Why do we find a dispersion that is significantly broader 
than Poisson even when 02 is only 2% larger than unity? Let 
us re- write a 2 as 



2 1 , a 

Q-2 = 1 + 



(AT) ^ (TV) 



1 -, 1 I a ' , 
1 + — I — -1 



(TV) 



(11) 



These studies, which have usually been done in the context 



From this expression, it is clear that any distribution where 
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a increases more slowly than (TV) will lead to Q2 ~ 1 for large 
values of (TV). Equation (111 also shows that the deviation 
from «2 = 1 can be directly related to the fractional scatter 
si in the model of equation dsj) : 



si = 



1 



(12) 



Even slight departures from «2 = 1 are therefore important: 
for example, if 02 = 1.1, then si = 0.46, i.e., a 46% intrinsic 
scatter. These results suggest that more care is needed in 
treating the scatter in HOD modeling. It is still true that 
aj 3 specifies the error incurred by assuming Poisson statis- 
tics when computing the j moment of the HOD, however. 
Using Poisson statistics will therefore result in ~ 5% error 
in the second moment and ~ 8% error in the third moment 
of the HOD. 

The subhalo occupation distribution can be well- 
approximated by the Negative Binomial distribution, which 
is given by 



P(N\ r, p) = 



r(r)r(jv + i) 



p r (i-p) 



(13) 



Here T(x) = (x — 1)! is the usual Gamma function and the 
parameters r and p are determined by the mean (TV) and 



variance a of the distribution: 

(TV) 2 



(TV) 

rr2 



(TV) 



(14) 



Using our model for the variance in TV, the values of p and r 
can be expressed in terms of (TV) and the fractional intrinsic 
scatter s\ alone: 



V = 



1 



1 



1 + s 2 (TV) 



(15) 



Given (TV), it is therefore straightforward to compute 
the full distribution P(TV| (TV), s 2 ). Note that in the limit 
r — ► 00 - i.e., as the intrinsic scatter goes to zero - the 
Negative Binomial distribution approaches the Poisson dis- 
tribution with mean (TV), as desired. 

The good agreement between the measured subhalo oc- 
cupation distribution and the Negative Binomial distribu- 
tion is illustrated in Fig. |10| Black data points show the sub- 
halo occupation distribution for (TV(> fj,)) = 4 (top panel, 
corresponding to p > 2.3 x 1CP 3 ), 15 (middle panel, corre- 
sponding to p > 5.6 x 10~ 4 ), and 125 (bottom panel, cor- 
responding to n > 7.1 x 10~ J ). The upper two panels use 
MS-II hosts with 10 12 < M vir < 10 12 ' 5 h' 1 M , while the 
lower panel uses MS-II hosts that are a factor of ten more 
massive. Poisson (green lines) provides a fairly good match 
for (TV) = 4 (the dispersion is 9% smaller than that of the 
data) but is noticeably too narrow at (TV) = 15 (the Pois- 
son dispersion is 25% smaller than that of the data). For 
(TV) = 125, Poisson is a very poor match to the data, with 
a dispersion that is too small by a factor of 2.25. The Neg- 
ative Binomial distribution (magenta curves) matches the 
data well for all three values of (TV). 

While the results presented in this section were derived 
using the cumulative mass function, we have confirmed that 
they hold equally well for the cumulative maximum circu- 
lar velocity function. Specifically, the intrinsic scatter in 
(TV(> Vmax/K/ir)) is nearly identical to that in (TV(> fx)}, 
i.e., the subhalo abundance also shows an intrinsic scatter 
of 18% when expressed in terms of Vma.x/V v i T . This is not 
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Figure 10. The subhalo occupation distribution from host ha- 
los in the MS-II with 10 12 < A/ vir < 10 12 5 h~ x M e (upper two 
panels) and with 10 13 < A/ vir < 10 13 5 h _1 Mq (bottom panel). 
Black data points show the measured distribution, green lines 
show the Poisson distribution with the same (TV), and magenta 
lines show the Negative Binomial distribution with the same (TV) 
and variance given by equation |9jl with fractional intrinsic scat- 
ter si = 0.18. Both distributions match the data for (TV) = 4 
(top panel, corresponding to [i > 2.3 X 10 — 3 ), while the Pois- 
son distribution is noticeably too narrow for (TV) = 15 (middle 
panel, corresponding to fi > 5.6 X 10~ 4 ) and much too narrow 
for (TV) = 125 (bottom panel, corresponding to /1 > 7.1 X 10 — 5 ). 
The Negative Binomial distribution matches the data well at all 
(TV). 
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true when considering Vmax/Vmax.host, however. This defini- 
tion introduces additional scatter in the subhalo abundance 
due to the scatter in concentrations at fixed mass: host halos 
of a given l/ v i r have a range of V max values. 



Ishiyama et al. ( 2009 1 have recently studied a sample of 



125 well-resolved halos from a cosmological simulation; their 
work corroborates this result. They show that the scatter in 
N(> Vmax/Vmax.host) is of order 40% for host halos with 
10 12 < A/vir < 2 x 10 12 h" 1 Mq. They find the scatter to 
be markedly reduced when using a criterion more similar 
to ours, N(> Vmax/V^oom.host) for subhalos within i?200m, 
signifying that subhalo properties correlate more tightly to 
host V v h than to host V ma x (see also Springel et al. 2008). 
This is yet another reminder that results on subhalo mass 
and circular velocity functions are sensitive both to how the 
host halo sample is defined and to how subhalos are selected 
(e.g., the limiting radius chosen). 

Finally, we have investigated whether the subhalo mass 
function shows systematic variation with the host's envi- 
ronment (as measured by the overdensity on a scale of 
5/i _1 Mpc) or host halo properties. No obvious correlations 
exist. Furthermore, the variation in the subhalo mass func- 
tion between hosts can be modeled as a simple normalization 
shift, i.e., the mass function within individual halos is con- 
sistent with the low-mass slope measured from the ensemble 
average. This is in good agreement with the behavior of the 
mass Aquarius halos' mass functions (e.g., Fig. [5J . 



4.3 Massive subhalos 

The parameter JH in equation |7|l is approximately the mass 
fraction relative to M v i r of the most massive subhalo in a 
"typical" host; /Ii = 0.01 therefore shows that the most mas- 
sive subhalo in a Milky Way-mass halo typically has about 
1% the mass of its host. The distribution of /j, for the most 
massive subhalo in each of our hosts is shown in Fig. |ll| Data 
are plotted for both (ii = M su b,o /Af v i r ,o (solid histogram and 
filled squares) and for /xi = Afa C c/M v i r ,o (dashed histogram 
and open squares^. Values for the individual Aquarius ha- 
los are shown by colored symbols. Fig. [TT] confirms that the 
probability distribution function for /xi peaks at /xi w 0.01 
for subhalo masses measured at z = 0, albeit with a large 
spread: the probability for m to be larger than 0.025 is 29% 
and to be larger than 0.1 is 4.3%, while 4.8% of halos have 
fii < 0.0025. 

The model presented in Sec. |4.2| suggests the scatter in 
the subhalo mass function should be nearly Poissonian for 
masses relevant for the most massive subhalo. We can test 
this directly: if the most massive subhalo in each halo is Pois- 
son sampled from the underlying (iV(> /x)), then the proba- 
bility density P(fJ,i) for finding /xi in the range [/xi, /xi+d/xi] 
can be directly computed as 



P(/xi) d/xi 



d/xi 



-<JV(>Ml)> 



dxti 



(16) 



11 



The cyan lines in Fig 
M au b,o (solid line) and M-, 



show P(m) computed this way for 
(dashed line); the agreement 



6 Note that what is plotted as the dashed histogram in Fig. |ll| 
is the probability distribution for the subhalo at z = with the 




Figure 11. Distribution over log^i of the most massive subhalo 
within -f?vh- of each of our Milky Way-mass hosts, both in terms of 
A^sub (solid histogram) and M aC c (dashed histogram). In each 
case, fii is computed relative to M v i r o- The values of fii for each 
Aquarius halo in the MS-II are also noted with square symbols 
(filled for Afgub^Oi open for M aC c). The distribution of fi± is fairly 
broad and peaks at /Xi « 0.01 (0.035) for M subj0 (M acc )- The 
cyan lines show the prediction for n\ (solid) and /J a cc,l (dashed) 
if the most massive subhalo is Poisson-sampled from our analytic 
fit to the cumulative mass function. In both cases, the Poisson 
prediction agrees very well with the actual distribution. 



with the actual distributions for /xi (the histograms) is ex- 
cellent. The distribution of zxi is therefore consistent with 
the Poisson hypothesis, both for the redshift zero masses 
and for the accretion masses. 

The distribution of xxi computed with respect to M acc 
can be used to estimate the probability that the Milky Way 
should host a galaxy at least as massive as the Large Magel- 
lanic Cloud (LMC), the most luminous MW satellite. Using 
the abundance matching results of |Guo et al. | (|2010[> and 
a ste llar mass of 2.5 x 10 9 A/© for the LMC (|Kim et al. 
1998| ), we find that M acc (LMC) 
Milky Way mass of 10 12 M< 



1.9 x 10 11 Mq. For a 
, there is a 8% chance of having 



a satellite with this accretion mass or greater. For a mass 
of 2.5 x 10 12 Mq - consistent with the abundance matching 
results of Guo et al. and the Local Group timing argument 



value obtained by Li & White (20081, this probability rises 
to 27%. 

This same line of reasoning can be used to compute 
the probability of a dark matter halo hosting two satel- 
lite galaxies at least as massive as the Small Magellanic 
Cloud (SMC). The SMC has a total stellar mass of ap- 
proximately 3 x 10 8 Mq ( Stanimirovic et al. 20041, corre- 



largcst M a 



not Af acc of the most massive subhalo at z = 0. 



sponding to an abundance matching mass at accretion of 
8.2 x 1O 1O M . If the Milky Way has a dark matter halo 
of mass 10 12 Mq, then there is a 3.3% chance of having a 
second-ranked satellite with at least this infall mass. For a 
halo mass of 2.5 x 10 12 Mq, the probability becomes 20%. 

These probabilities ignore any scatter in the M t — Afhaio 
relation. This scatter is subdominant to the uncertainty in 
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the Milky Way's mass, however. For example, |Guo et al.| 
(20101 give 80% confidence intervals that correspond to 



±25% on the infall halo masses of the Magellanic Clouds 
based on a dispersion of 0.2 in Iog 10 (Mhak>) at fixed M*. 
Clearly, the satellite statistics quoted above give some sup- 
port for the higher estimates of the Milky Way's halo mass. 

4.4 Subhalo accretion times 

Subhalos persisting to z = are fossils of previous halo 
merger events. The accretion redshifts z acc of z = sub- 
halos are therefore important quantities for understanding 
the dynamical evolution of subhalos and can help constrain 
other quantities, such as merging timescales. The median 
z acc for z — subhalos as a function of /j, = M su b,o/M v i r is 
shown as a solid curve in Fig. |12| The 25th and 75th (10th 
and 90th) percentiles of the z acc distribution are included 
as dashed (dotted) curves. Results for Milky Way-mass host 
halos are shown in the left panel. The median z acc is a de- 
creasing function of [i, which is mainly a dynamical effect: 
massive subhalos that are accreted at early times are able to 
lose their orbital energy and angular momentum and merge 
with the dominant subhalo by z = while low-mass subha- 



los are not (e.g., Boylan-Kolchin et al. 20081. As a result 



the only remaining massive subhalos are those that were ac- 
creted at relatively late times. At fi = 10~ 3 , the typical z acc 
of a surviving subhalo is somewhat larger than z = 1, while 
at fi = 10 _1 , Zacc ~ 0.3 for surviving subhalos. It is diffi- 
cult to compare the individual Aquarius halos, as there are 
only 60 satellites across all six halos satisfying /i > 0.001 
and d < R v i r , but the median accretion time for these 60 
satellites, z acc = 0.76, falls well within the expected range 
based on the distribution from the MS-II. 

For comparison, the right panel of Fig.[l2]shows the dis- 
tribution of z acc for surviving subhalos in much more mas- 
sive hosts, those with M vir > 10 13,5 h' 1 M®. The trend is 
very similar to that seen for Milky Way-mass hosts, namely, 
the median z acc is a decreasing function of /i. The typical 
z a cc values are uniformly lower for the massive halo sample, 
however: at fi = 10 -3 , the median surviving subhalo in the 
massive host sample was accreted at z acc « 0.6. This dif- 
ference is a manifestation of the hierarchical growth of dark 
matter halos: more massive halos are dynamically younger 
and have accreted a larger fraction of their mass (and of their 
subhalo population) recently. Fig. 



12 



also indicates that a 
large fraction of the galaxy population in galaxy clusters has 
joined the cluster fairly recently, whereas most of the lumi- 
nous satellites of galaxies similar to the Milky Way have 
been part of their host halos since z « 1. At all masses, the 
median value of M a . 
that a typical z 



c/M su b,o is 2.4 (see Fig. All, indicating 
: subhalo has lost 60% of its mass. 
The results in the left panel of Fig. |12| app ear to differ 
substantially from those of Gao et al. (20041, who found 



that 90% of subhalos in a Milky Way-mass halo at z = 
have z acc < 1. These differences are due to the choice of 
definition for z acc . Gao et al. define z acc as the last time a 
subhalo entered the FOF group of the main progenitor of its 
z = host while we define z acc to be the redshift when the 
subhalo's mass was at its maximum. This latter definition 
is closer to the first time a subhalo enters its host's FOF 
group. We have checked that using Gao et al.'s definition 
of z acc gives results similar to theirs. Our results for more 



massive halos (the right panel of Fig. 12 1 agree with the 
findings of De Lucia et al. (2004 e.g., their fig. 10), who use 
the same definition of z acc as we do. 



5 MERGER HISTORIES 

Galaxy-mass dark matter halos are constantly bombarded 
by smaller halos in the ACDM model. This has been the 
source of great concern with respect to the formation of 
MW-like galaxies: how is it possible to reconcile an active 
merging history with the existence of a thin stellar disk (e.g., 
Toth fc Ostriker| [1992] |Velazquez 



White||1999| |Benson| 
et al.|2004||Kazantzidis et al.|2008||Hopkins et al.|2009| )? In 

order to understand the severity of the problem, we investi- 
gate the merger histories of MW-mass halos in this section, 
concentrating on mergers of objects with mass comparable 
to that of the MW disk. 

When considering the merger history of a dark matter 
halo, two major questions are (1) what is the rate of infall 
of other, smaller dark matter halos across the halo's virial 
radius and (2) what is the rate at which other dark matter 
halos lose their identity and merge with the center of the 
halo? Clearly, the set of objects that merge into the center 
of a halo will be a subset of those that crossed the halo's 
virial radius at earlier times, but connecting the two is far 
from trivial. Several recent papers have addressed the fre- 



quency of halo- halo mergers. Fakhouri & Ma (20081, Guo 



fc White| ( |2008| ), and |Genel et al.| ( |2009| ) investigated the 
halo merger rate directly in the MS, while [Neistein fc Dekel| 
(20081) andlCole et al.|(|2008|) used extended Press-Schechter 



theory to derive merger rates and compared the results to 
the MS. The most relevant study for this paper is |Stewart| 
( 2008 1 , who used an iV-body simulation to extract 



et al. 



halo-halo merger rates for MW-mass halos. Stewart et al. 
found that while recent major mergers are very rare, merg- 
ers involving halos with masses of ~ O.lM^o are quite 
common over the past 10 Gyr and mergers with halos hav- 
ing masses exceeding the mass of the galactic disk are vir- 
tually inevitable over the same time period. These results 
led the authors to conclude that the frequency of mergers 
may present a substantial challenge to our understanding 
of disk galaxy stability, as the majority of galaxies in MW- 
mass halos are disk-dominated ( Weinmann et aL||2006 van 



den Bosch et al.|2007[|Park et al.|2007| ) 

The mass resolution of the MS-II allows us to expand 
on the results of Stewart et al. by studying mergers with the 
central regions of MW-mass halos. These are more clearly 
relevant to the survival of stellar disks (which have charac- 
teristic radii that are less than 10% of -Rvir). Accordingly, 
in this section we consider only mergers between accreted 
subhalos and the dominant subhalos in MW-mass halos. We 
construct samples of merging objects by searching the main 
progenitor branches of the dominant subhalos for merger 
events with subhalos satisfying V aC c > 50kms _1 . 

In order to ensure that the mergers we are considering 
are well resolved and are not affected by limited numeri- 
cal resolution, we have computed the orbits and internal 
properties of subhalos immediately prior to merging. For the 
mergers considered here - those with V acc > 50kms _1 for 
the satellite - the median separation between the satellite 
and the dominant subhalo immediately prior to merger is 
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Figure 12. Rcdshift of accretion 2 acc as a function of M su b,o/Af v j r for all resolved subhalos in galaxy-mass (left) and cluster-mass 
(right) host halos. The median relation is shown as a solid curve while the 25 and 75 (10 and 90) percentiles are shown as dashed (dotted) 
curves. More massive subhalos at the present day were accreted later on average than low-mass subhalos; the spread of accretion times for 
massive objects is also smaller. Subhalos of more massive host halos are accreted at lower redshifts on average, reflecting the hierarchical 
build-up of dark matter halos. 



about 15ft 1 kpc; in 80% (90%) of the cases, the separation 
is less than 30 (50) h~ kpc. 

Furthermore, the satellites themselves are well-resolved 
immediately prior to merging: as we show in the Appendix, 
even halos at our lower limit of 14cc = 50 km s _1 have on av- 
erage 70 particles just prior to merging, while satellites with 
Vacc = 65kms _1 typically have 150 particles at the time 
of merging. Accreted halos that are capable of substantially 
impacting the Galactic disk can therefore be reliably tracked 
to within < 10% of i? v ir, allowing us to separate true merg- 
ers with the centers of halos from accretion events that cross 
i?vir but remain at large halo-centric distances. 

Figure [13] explores the probability for dominant 
z — MW-mass subhalos to have merged since 
redshift z with another subhalo with M aC c/Af v ir,o > 
0.01,0.03,0.05,0.07,0.1, and0.15 (top to bottom). Approx- 
imately two-thirds of MW-mass halos have experienced a 
merger with a halo having M acc /M v j r: o > 0.03 since 2 = 1, 
and over 90% have had at least one such merger since 2 = 3. 
On the other hand, only 30% of halos have merged with a 
■Macc/A^vir.o > 0.1 halo since 2 = 1 and just 50% of halos 
have ever experienced such a merger. Approximately 50% 
of halos have had a merger with 0.03 < M acc /M v i ri o < 0.15 
since z=l. 

It is important to note that we have been quoting satel- 
lite masses in terms of the maximum dark matter mass the 
satellite has ever attained, M acc - Subhalos usually lose a 
large fraction of their mass before merging: for our sample, 
the average subhalo loses 90% of its mass between accretion 
and merging (see Fig. |A2[ ), independent of M acc . Subhalos 
with M^cc/Mvirfi = 0.1 therefore typically have a mass ratio 
of 1:100 at merger, meaning the mass involved in a potential 
disk impact may be three times smaller than the disk mass 
rather than three times larger. 

We can also investigate the distributions of merger his- 
tories in terms of the n th -most massive halo with which a 
MW-mass halo has merged. This is done in Fig. |14| The 
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Figure 13. The fraction of central subhalos that have merged 
with an object having Af acc greater than 1, 3, 5, 7, 10, and 15% of 
A^vir.o (upper through lower curves) since redshift z, as a function 
of z. Virtually all halos have merged with a M acc > 0.01 M v i T q 
halo since 2 = 3, while only half have merged with a M acc > 
0.10 M v ; r o halo in the same redshift range. Recent merger events 
with M acc > 0.03 Af v i r o are common: 70% of halos have had such 
a merger since 2 = 1. 

solid lines show the probability of having a merger with 
Af a cc/Afvir, o > M as a function of n, while the dashed curves 
give the same probability for mergers at z < 2. Different 
colors correspond to different merger ranks: green shows the 
probability distribution for the most massive merger for each 
halo (in terms of A/ acc ) while red shows the distribution for 
the 5th most massive merger. 90% of halos have had one 
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Figure 14. The cumulative probability of having a merger with 
fi > M acc /M v i rj o for the 1st, 2nd, 3rd, 4th, and 5th-most massive 
central merger that a given halo experiences, both over the halo's 
entire lifetime (solid lines) and since z = 2 (dashed lines). The 
median value for the most massive merger is fj, = 0.1, while the 
3rd-most massive merger typically has fj, = 0.02. 



merger event with M acc > 0.03M v i rj o since z — 2, 60% have 
had two such events, and 30% have had three. The median 
value for the most massive merger a halo has experienced 
since z = 2 is M acc /M v i ri o = 0.1, and the median redshift for 
such merger^] is z = 0.75; only 10% of halos have had two 
events with M acc /M v i ri o > 0.1. The most massive merger 
event since 2 = 2 for each of the Aquarius halos (A-F) has 
Af acc /M vir , of 0.033, 0.078, 0.028, 0.030, 0.033, and 0.259, 
respectively. Halos A, C, D, and E therefore lie near the 10% 
most quiescent halos of Fig. [14] while halo B is fairly typical 
and halo F has had a recent major merger, which is rare. 

A thorough understanding of the stability of galactic 
disks in the presence of satellite infall requires the inclu- 
sion of baryonic physics, as several studies have shown that 
factors such as the orbits and gas content of the merging 
galaxies can have a substantial impact on remnant prop- 
erties JSpringel & Hernquist 2005 ; |Kazantzidis et aT||2008| 



2009; Scannapieco et al.||2009| |Hopkins et al. 
et al.||2009| |Stewart et al.||2009| |Moster et a 


|2009| |Purceli 
J. |2010| |Gov- 


ernato et al. 


2009 


D'Onghia et al. 


2010 


) Nevertheless, we 



can still use the statistics computed above to learn about 
the frequency of satellite-disk interactions, as the dynamics 
that determine which subhalos merge are mostly set by dark 
matter halo properties such as the satellite's orbit and mass 
at infall. 

Approximately 90% of MW-mass halos have had a cen- 
tral merger with a halo of M acc > 0.03 M v i r ,o since z = 2. If 
such mergers destroy galactic disks, it is extremely difficult 
to understand the frequency of disk galaxies in MW-mass 

7 This is the redshift at which the satellite subhalo merges with 
the central subhalo, not the redshift at which the progenitor halos 
merge. 
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halos (g 70%; |Weinmann et al.||2006 Ivan den Bosch et al 



2007| |Park et al.|2007[ ). Around half of MW-mass halos have 
experienced a merger of M acc > 0.1 M v i ri o since z = 2. These 
events are likely to strongly affect galaxy disks; however, the 
frequency with which they occur, coupled with the possibil- 
ity that a non-negligible fraction of such mergers may have 
sufficient gas to reform a disk in the merger remnant, does 
not seem high enough to present a major obstacle to forming 
galactic disks in a majority of MW-mass halos. It is there- 
fore essential to simulate a large number of mergers in the 
mass range of 0.03 < M acc /M v i r ,o < 0.1 with a distribution 
of realistic orbits and galaxy models, including a variety of 
gas fractions. This will significantly advance our knowledge 
of disk formation and survivability in the ACDM cosmology. 



6 DISCUSSION 

6.1 How (a)typical are the Aquarius halos? 

In general, the agreement between the properties of the 
Aquarius halos and those of the full sample of MW-mass 
halos from the MS-II is quite good: no halo is an outlier 
from all of the relations studied here. Halos A and C form 
somewhat earlier than average. Halo C also has a larger con- 
centration than is typical. On the other hand, neither A nor 
C seems to be an outlier in the distribution of the spin pa- 
rameter A'. Halo F is somewhat unusual in that it has a 
recent major merger; on the other hand, it is quite typical 
in most of its other properties. The Aquarius halos seem to 
reflect the diverse properties of Milky Way-mass halos, both 
in terms of assembly history and 2 = structure. 

The Aquarius halos do not, however, seem to uniformly 
sample the spin parameter A': five of the six halos lie below 
the median of the distribution, while the sixth is a merger 
remnant and may well relax to a lower value as well. We 
have explicitly checked that the isolation criterion in the 
Aquarius halo selection does not bias the spin parameters, 
as the distribution of A' from the isolated and non-isolated 
samples are statistically identical. The somewhat low spin 
parameter distribution of the Aquarius halos appears to be 
a statistical fluctuation. This may be connected to the rel- 
atively quiet recent merger history of the central subhalos 
for the Aquarius halos. 



6.2 How (a)typical is the Galaxy's halo? 

While the Milky Way is frequently considered as the proto- 
typical massive spiral galaxy, this has its roots as much in 
convenience as in evidence. While it is difficult to use dissi- 
pationless iV-body simulations to discern whether the MW 
is indeed a typical galaxy for a halo of its mass, we can in- 
vestigate whether the Galaxy's halo itself is typical among 
those of similar mass. Many possible constraints are pred- 
icated on a more precise determination of the mass of the 
MW's halo, however. Current estimates vary by a factor of 
approximately 2 to 3 (with 1 < M v j r ,MW < 3 x 1O 12 M0). 
This uncertainty has a large effect on the probability of host- 
ing massive satellite galaxies such as the SMC and LMC: the 
chances of having a subhalo capable of hosting the LMC or 
two subhalos capable of hosting the LMC and SMC are ap- 
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proximately 3-8% for a MW mass of 10 12 Mq but rise to 
20-27% if the MW's halo has a mass of 2.5 x 10 12 M . 

On the other hand, a more massive Galactic halo implies 
that merger events of a given M acc /M v j rj o correspond to a 
larger fraction of the MW's disk mass. If the MW's thin disk 
must be explained via an unusually quiescent merger history 
over the past 10 Gyr (e.g., Hammer et al. 20071, we find 



tude of JV (> Knax/Knax.host) between a s = 0.74 and 1.0 is 
only 14% dReed et al 120051 iDiemand et al.||2008| iDiemand 



that the MW's halo must be even more unusual (in terms of 
merger history) if the halo is massive. If we assume the MW 
cannot have merged with a M acc > 3Afdi s k,MW halo since 
z = 2, this corresponds to fi = 0.1 for Mmw = 1 12 Mq and 

12 Mr, 
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we find 



/i = 0.04 if M MW = 2.5 x 10 12 Mq. From Figure 
a 50% chance for the former case but only a 20% chance 
for the latter. If halos with M acc = Mdisk.MW destroy thin 
disks, the probabilities drop to 12% and 4%, respectively. 

A possible resolution of this tension - that the existence 
of the LMC and SMC argue for a more massive MW halo 
while the lack of recent mergers argues for a less massive 
Galactic halo - is to directly connect the two issues. A re- 
cent analysis of the LMC's orbital history ( Besla et al.|2007 l 
based on the updated proper motions of the LMC ( |Kallivay-| 
alii et al. 2006) indicates that the LMC likely fell into the 
Milky Way's halo within the last ~ 3 Gyr. The LMC's in- 
fall mass is likely to be in excess of 10 1 h~ Mq (and the 
SMC's infall mass is likely to be only slightly less massive; 
Guo et al.|2010 1, meaning that the Milky Way has actually 
experienced a recent 1:10 merger at the halo level. [Stewart] 
et al. ( 2008 I show that Milky Way-mass halos typically un- 



dergo one such merger at the halo level over the past 10 Gyr, 
perhaps indicating that the Milky Way's accretion history 
is special only insofar as the LMC fell in so recently. 

6.3 Cosmological Parameter Dependence 

As noted in Sec. |2.1| the cosmological parameters used for 
the MS-II are slightly different from the most recent obser- 
vational estimates: the values of erg = 0.9 and n s — 1.0 are 
approximately 3<r away from the best-fitting values of 0.81 
and 0.96, respectively. Many of our results are insensitive 
to such changes. The abundance of MW-mass halos differs 
by only ~ 7% in the WMAP5 cosmology compared to the 
Millennium-II cosmology, for example, and the basic prop- 
erties of halo assembly histories scale weakly with erg 



den Bosch|2002[ ). The subhalo abundance of MW -mass halos 
could potentially be affected by as and n 3 , as decreasing as 
tends to reduce halo and subhalo concentrations at a fixed 
mass (e.g., |Maccio et al.|2 008). This would tend to decrease 
the amplitude of {N(> fi)} in lower as cosmologies, as sub- 
halos would be more easily disrupted. On the other hand, 
halos tend to form later in such cosmologies, meaning that 
subhalos are subject to less dynamical evolution (see the dis- 
cussion at the start of Sec. 4.1 1; this would tend to increase 



the amplitude of {N(> fi)) for a lower erg. These two effects 
compensate to some extent, reducing the overall shift. 

In order to properly ascertain the effects of cosmology 
on the structural properties of MW-mass halos, it is neces- 
sary to run simulations that differ only in their cosmological 
parameters and to search for differences in statistical prop- 
erties. Nevertheless, there is reason to think that changing 
from the cosmology of the Millennium, Millennium-II, and 
Aquarius simulations to the WMAP 5 cosmology will not 
have a strong impact. For example, the difference in ampli- 



|fc Moore||2009[ ), so reducing a 8 from 0.9 to 0.81 will likely 
affect the amplitude of the subhalo mass function by a few 
percent only. 



7 CONCLUSIONS 

We have presented a statistical study of Milky Way-mass ha- 
los from the Millennium-II Simulation, which contains over 
7000 halos in the mass range 10 11 ' 5 < A/ vir < 10 12 ' 5 h" 1 Mq. 
Our principal results can be summarized as follows: 

• As several previous studies have shown, halo growth 
proceeds in an "inside-out" fashion: the central gravitational 
potential of a typical MW-mass halo reached half of its 
present-day value by z — 4, on average, whereas the virial 
mass reached half of its present-day value at z — 1.2. 

• The ratio of z — mass to the host halo mass for the 
most massive subhalo of a MW-mass halo has a broad distri- 
bution that peaks at 1%. The corresponding quantity com- 
puted using infall mass M acc instead of M su b,o peaks at 3.5% 

• The differential and cumulative subhalo mass functions, 

computed in terms of fi = M su b/M v ir, are both well-fitted 

by a power law with an exponential cut-off at large /x. We 

find both the abundance per log decade in /i and the number 

of subhalos with mass greater than [i to be proportional to 
^-0.935 £ Qr sma rj ^ 

• The scatter in the cumulative mass function at fixed fi 
is only well-approximated by a Poisson distribution at large 

(corresponding to mean occupation number (TV) < 4 — 5). 
Intrinsic scatter of approximately 18%, independent of the 
host halo mass, becomes increasingly important at lower fi 
(larger (TV)) and is dominant over Poisson scatter at [i < 
3 x 10~ 4 ((N) > 20). 

• The Negative Binomial distribution with variance in 
{N(> fi)} that is equal to the sum of a Poisson term and a 
fractional intrinsic scatter of 18% matches the subhalo occu- 
pation data well at all /i for host halo masses between 10 12 
and 10 15 h,- 1 M G 



z . 



• The statistic Q2 (m) [equation ( 10 1] , which is frequently 
used to characterize deviations of the HOD from Poisson, 
is not discriminating: distributions that differ strongly from 
Poisson can lead to few percent differences of a.2 from unity. 

• Accretion redshifts of z = subhalos vary systemat- 
ically with host halo mass and with the ratio of subhalo 
to host masses. Massive subhalos at z = were typically 
accreted more recently. This is a dynamical effect: massive 
subhalos accreted at early times either merge via dynamical 
friction or are heavily tidally truncated. At fixed fi, subhalos 
of more massive host halos were accreted more recently, re- 
flecting the hierarchical build-up of dark matter halos. Sub- 
halos with fj, = 0.001 have a median accretion redshift of 
Zucc = 1-1 for hosts with M V1T ~ 10 12 ft -1 Mq and z acc = 0.6 
for halos with M vir > 10 13 5 Mq. 

• The frequency of central mergers - that is, the merger 
of a satellite with the central region of a MW-mass halo 
- is a strong function of the satellite's M acc . Since z = 2 
90% of MW-mass halos have experienced a central merger 
with a satellite having M acc /M v i r ,o > 0.03 while 50% have 
had such a merger with a M acc /M v i r ,o > 0.1 satellite, and 
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only 30% with a M aC c/Afvir,o > 0.15 satellite. The compat- 
ibility of thin stellar disks with the frequent mergers ex- 
pected in ACDM thus depends strongly on exactly which 
mergers destroy disks. If M acc /M v ir,o > 0.03 events are suf- 
ficient, the Milky Way must have an accretion history that 
lies among the quietest 10% for halos of similar mass, while 
if M acc /Mvi r .o > 0.1 is required, then the Milky Way may 
have had a fairly typical merger history. 
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APPENDIX A: SUBHALO MASSES AND 
CIRCULAR VELOCITIES 

Figure |A1| shows the median relations between M acc and 
Vmax.o (solid black line) and between M su b,o and V m ax,o 
(solid blue line); the latter is identical to the one plot- 
ted in Fig. [I] Also plotted are the regions centered on the 
median M acc — V max ,o relation containing 50 and 80% of 
the distribution (dashed and dotted curves, respectively). 
The M acc — Vmax.o relation follows the same slope as the 
-M au b,o ~ Vrnax.o relation but has an amplitude that is ap- 
proximately 2.4 times larger: the solid magenta line shows 
the result of multiplying the fit to M su b,o — Kn ax ,o at z = 
by a factor of 2.44. The average subhalo at z = was there- 
fore approximately 2.4 times more massive at accretion, irre- 
spective of its present V m ax- The spread in the M acc — V^ax.o 
relation is non-negligible: the 90% value of M acc is typically 
a factor of 2.5 larger than the 10% value at fixed Kn ax . 10% 
of halos have lost over 80% of their mass while another 10% 
have lost less than 30%. 

Figure [A"2] shows the median relation between M acc and 
V a cc (solid black curve), along with the ranges containing 
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Figure Al. M acc as a function of V m ax for subhalos surviving 
to z = 0. The solid curve shows the median M acc , while the 
dashed (dotted) curves show the 25 and 75 (10 and 90) percentiles. 
Also plotted, in magenta, is the best-fitting power law relation 
between M acc and V max ,o, as well as the median relation between 
M(z = 0) and V max ,o for subhalos (blue curve). The virtually 
constant offset between the blue and solid black curves shows 
that the typical surviving subhalo loses a factor of ~ 2.5 in mass 
between infall and the present day, irrespective of Vin a x,o- 



Figure A2. Af acc as a function of Vm ax for subhalos with 
V aC c > 50kms that merge with their host. The solid curve 
shows the median M acc , while the dashed (dotted) curves show 
the 25 and 75 (10 and 90) percentiles. Also plotted, in magenta, is 
the relation between Vm ax and M merge , the mass of the subhalo 
immediately prior to merging. The median V max — M merge rela- 
tion is nearly parallel to the median V max — M acc relation but is 
lower in amplitude by approximately 1 dex, indicating that sub- 
halos lose approximately 90% of their mass before merging for all 



50 and 80% of the distribution (dashed and dotted black 
curves). Our lower limit of Va CC = SOkms -1 corresponds to 



a median M acc of 5 x 10 9 ft -1 M Q . Figure A2 also shows in 
magenta the relation between Vacc and M morge , defined to be 
the mass of the subhalo at the snapshot immediately prior 
to merging. The relation is nearly parallel to M acc (V acc ) (the 
solid black line) but is lower in amplitude by approximately 
a factor of 10 at all V acc - This means that these subhalos 
lose approximately 90% of their mass prior to merging, ir- 
respective of their mass at infall. 
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